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INTRODUCTION 


The  flow  and  heat  transfer  in  the  projectile  launching  tube  of  a  weapon 
is  typically  a  complicated  two-phase  flow  where  combustion  products  are  mixed 
with  unburned  propellant  grains.  A  detailed  calculation  of  the  flow  field  in 
the  gun  tube  would  provide  important  information  such  as  local  transient  heat 
transfer  rates  and  propellant  burning  characteristics.  This  information  would 
contribute  to  the  understanding  and  solution  of  problems  associated  with  gun 
barrel  erosion  and  catastrophic  gun  failures.  The  present  work  is  an  extension 
in  the  development  of  the  computer  code  ALPHA  (Ref.  1)  which  solves  the  governing 
equations  for  the  two-phase,  two-dimensional  flow  in  a  gun  tube. 

Other  efforts  in  modeling  the  flow  phenomena  in  guns  include  quasi-one- 
dimensional  inviscid  two-phase  flow  analyses  of  the  propellant  combustion  process 
(e.g.,  Ref.  2-7),  a  two-dimensional  inviscid  two-phase  flow  analysis  (Ref.  8), 
and  time-dependent  boundary  layer  analyses  applied  to  the  flow  of  propellant 
gases  in  a  gun  tube  (Refs.  9  and  10).  The  boundary  layer  procedures  suffer  from 
the  shortcoming  that  the  starting  conditions  near  the  projectile  base  are  not 
well  defined,  and  according  to  conventional  boundary  layer  theory  the  heat  flux 
near  the  base  approaches  infinity  because  the  base  is  a  singular  point  (Ref.  10). 
Furthermore,  the  validity  of  the  boundary  layer  approximations  is  questionable 
at  both  the  breech  end  and  the  projectile  base  region,  and  even  the  most 
sophisticated  boundary  layer  analysis  presently  used  for  gun  barrel  problems, 
e.g.,  Refs.  9  and  10,  did  not  consider  the  two-phase  flow  aspects  of  the 
propellant  combustion  process.  The  significant  features  of  the  two-phase  flow 
interior  ballistics  codes  (Refs.  2-7)  were  reviewed  by  Kuo  (Ref.  11).  The  main 
objection  to  these  analyses  (Refs.  2-7)  would  seem  to  be  the  presumption  of 
quasi-cne-dimensional  flow  and  the  attempt  to  predict  heat  transfer  to  the  barrel 
using  rather  simple  unsteady  boundary  layer  models  or  correlation  formulas. 

The  two-dimensional  analysis  of  Gough  (Ref.  8)  utilizes  an  explicit  two-step 
time-dependent  algorithm  due  to  MacCormack  (Ref.  12)  with  upwind  spatial  differ¬ 
encing  modifications  suggested  by  Moretti  (Ref.  13).  The  coordinate  transfor¬ 
mation  technique  used  in  Ref.  8  is  based  on  the  method  of  Thompson,  et.  al. 

(Ref.  14).  The  initial  work  of  Ref.  8  considers  only  a  stationary  coordinate 
system  to  model  the  ignition  process  in  a  closed  chamber  with  geometrical 
complicat ions  such  as  a  hemispherical  breech  closure  plug  and  a  boattail 
projectile  intrusion.  According  to  Ref.  8,  this  stationary  coordinate  system 
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would  be  transformed  to  a  moving  system  via  a  solid  phase  Lagrangian  transfor¬ 
mation  so  that  the  mesh  moves  with  the  particles.  It  is  not  known  if  this 
technique  will  give  accurate  results  when  projectile  motion  is  included  in  the 
calculation.  However,  since  the  governing  equations  in  Ref.  8  were  written  in 
nonconservation  form,  significant  coordinate  motion  may  give  rise  to  mass 
sources  due  to  time  truncation  error  introduced  by  a  moving  mesh  as  demonstrated 
by  Thomas  and  Lombard  (Ref.  15).  Furthermore,  if  the  mesh  motion  is  governed 
by  the  solid  particle  motion, inadequate  resolution  may  result  in  certain  areas 
of  the  computational  domain.  Finally,  prediction  of  barrel  heat  transfer  with 
an  inviscid  core  flow  analysis  would  require  potentially  inaccurate  unsteady 
boundary  layer  models  as  in  the  quasi-one-dimensional  analyses. 

The  complex  nature  of  the  flow  in  the  projectile  base  region  and  in  the 
breech  end  of  the  barrel  does  not  permit  simplifying  approximations  to  be  made 
in  the  governing  fluid  flow  equations,  and  therefore  in  principle  the  solution 
of  the  full  Navier-Stokes  equations  is  required,  rather  than  some  simpler 
approximate  set  of  equations.  Fortunately,  recent  developments  in  computational 
fluid  dynamics  have  made  possible  the  prediction  of  the  detailed  flow  field  in 
configurations  such  as  a  gun  barrel  using  the  full  Navier-Stokes  equations. 

The  equations  and  coordinate  system  developed  under  this  effort  have  been 
incorporated  into  an  existing  three-dimensional  time— dependent  compressible 
Navier-Stokes  calculation  procedure  (the  HINT  code)  which  was  originally 
developed  under  United  States  Navy  and  Air  Force  sponsorship  for  other  purposes 
by  staff  members  of  Scientific  Research  Associates,  Inc.  (Refs.  16-19).  The 
MINT  procedure  solves  the  governing  equations  using  a  consistently-split, 
linearized,  block-implicit  numerical  scheme  (Ref.  20). 

Under  a  previous  effort  (Ref.  1)  a  mathematical  model  of  a  two— phase, 
two-dimensional  flow  was  developed  and  a  computer  code  (ALPHA)  was  constructed 
for  the  numerical  solution  of  the  equations  resulting  from  this  mathematical 
model.  The  model  developed  consists  of  the  governing  equations  for  an  axisymmetr ic , 
two-phase  flow  in  a  gun  tube  with  a  rotating  projectile,  and  a  system  of  consti¬ 
tutive  relations  describing  the  molecular  viscosity  and  thermal  conductivity, 
turbulence  length  scale,  gas  equation  of  state,  intergranular  stress,  interphase 
drag,  interphase  heat  transfer,  and  solid  phase  combustion.  The  governing 
equations  and  corresponding  initial  and  boundary  conditions  described  the  firing 
cycle  beginning  with  a  fluidized  and  ignited  solid  phase,  and  ending  with  the 
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projectile  exiting  the  gun  tube.  Chemical  reactions  within  the  gas  phase  were 
excluded  from  the  formulation.  An  axisymmetric  time-dependent  adaptive  coordinate 
system  for  interior  ballistics  flow  field  calculations  was  developed,  and  the 
projectile  and  distinct  filler  elements  were  treated  using  a  quasi-one-dimensional 
lumped  parameter  analysis. 

Under  the  present  effort  the  ALPHA  code  has  been  modified  to  include  both 
an  ignition  model  to  permit  computation  of  the  complete  firing  cycle  and  the 
capability  to  treat  geometries  in  which  the  tube  radius  is  a  function  of  axial 
distance.  An  efficient  procedure  for  mass  storage  and/or  large  core  memory 
utilization  for  the  ALPHA  code  dependent  variable  array  and  the  block-tridiagonal 
matrix  inversion  array  was  implemented  and  made  operational.  Further,  an 
existing  SRA  capability  for  velocity  vector  and  contour  plotting  of  ALPHA  code 
output  was  made  operational  on  the  BRL  computer  system.  The  moving  coordinate 
system  capability  in  the  previously  developed  ALPHA  code  (Ref.  1)  was  observed 
to  have  truncation  error  produced  mass  sources  associated  with  the  form  of  the 
governing  equations  utilized.  This  observation  of  truncation  error  produced  by 
moving  meshes  is  consistent  with  the  findings  of  Thomas  and* Lombard  (Ref.  15), 
who  suggested  reformulation  of  the  governing  equations  in  strong  conservation 
form  along  with  a  special  technique  for  determination  of  the  local  time-dependent 
Jacobian  determinant  of  the  coordinate  transformation.  Therefore,  under  the 
present  effort  the  governing  conservation  equations  were  reformulated  in  a  strong 
conservation  form  by  application  of  a  Jacobian  transformation  to  the  equations 
in  cylindrical-polar  coordinates.  The  resulting  computer  code  has  been  desig¬ 
nated  as  ALPHA2. 
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THEORETICAL  ANALYSIS 


Approach 

The  governing  equations  for  a  two-phase  two-dimensional  flow  in  a  gun  tube 
which  were  originally  presented  in  vector  form  in  Ref.  1  are  repeated  below  for 
completeness.  The  provision  for  a  rotating  projectile  is  considered  by  solving 
the  azimuthal  momentum  conservation  equation  with  the  appropriate  boundary 
conditions  at  the  projectile  base.  The  governing  equations  may  be  obtained  by 
employing  either  the  time- averaging  procedure  utilized  by  Ishii  (Ref.  21)  or 
the  formal  averaging  approach  used  by  Gough  (e.g.,  Refs.  22,23)  or  Gough  and 
Zwarts  (Ref.  24).  In  the  present  derivation,  the  averaging  procedure  of  Gough 
(Ref.  22)  has  been  selected  because  of  its  notational  convenience;  however, 
extensive  reference  to  the  work  of  Ishii  (Ref.  21)  has  been  made  in  order  to 
verify  the  results  obtained.  In  the  following  analysis,  a  gas-solid  mixture 
is  assumed  with  a  constant  solid  phase  density,  p^.  Numerous  assumptions  and 
approximations  are  required  in  order  to  formulate  a  tractable  problem.  Most  of 
the  required  assumptions  have  been  stated  previously  by  Gough  (e.g.,  Ref.  2), 
and  those  necessary  in  the  present  work  are: 

(1)  The  gas  and  solid  phases  occupy  separate  complementary  regions,  and  within 
each  region  the  material  may  be  treated  as  a  homogeneous  continuum. 

(2)  The  flow  of  the  heterogeneous  mixture,  composed  of  the  two  interacting 
continua,  can  be  described  by  appropriately  defined  averages  of  the  flow 
properties . 

(3)  If  solid  phase  combustions  occurs,  the  energy  deposition  is  taken  to  be 
in  the  gas  only. 

(4)  The  solid  phase  is  deformable  and  incompressible.  However,  locally  no 
relative  motion  between  the  solid  particles  is  considered.  Thus  the  average 
stress  in  the  solid  phase  is  an  isotropic  normal  stress. 

(5)  The  influence  of  solid  phase  deformation  on  the  particle  surface  area  is 
neglected,  and  the  interfacial  average  of  the  particle  velocity  is  equal 
to  the  volume  average  in  the  absence  of  burning. 

(6)  The  interphase  drag  is  determined  from  steady  state  correlations;  the 
unsteady  virtual  mass  effect  is  not  considered. 

(7)  The  interphase  heat  transfer  is  determined  from  steady  state  correlations. 
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(8) 


The  Noble-Abel  equation  of  state  is  employed.  The  specific  heats 
(0^  and  c^)  are  taken  to  be  independent  of  temperature. 

(9)  The  regression  rate  of  the  surface  of  the  burning  propellant  is  a  function 
of  the  average  gas  properties  and  the  propellant  surface  temperature. 

(10)  Heat  transfer  to  the  solid  phase  is  treated  as  a  one-dimensional  process 
in  order  to  determine  the  propellant  surface  temperature. 

(11)  The  pressure  drop  at  the  gas-solid  interface  is  negligible. 

Governing  Equations 

Both  Ishii  (Ref.  21)  and  Gough  (Refs.  22,  23)  have  presented  the  relations 
for  the  average  of  time  and  space  derivatives  in  a  two-phase  mixture.  Using  the 
above  assumptions  a  system  of  partial  differential  equations  is  obtained  con¬ 
taining  interface-averaged  source  terms  arising  from  averaging  the  basic  conser¬ 
vation  equations  for  the  two-phase  mixture.  A  basic  quantity  used  to  describe 
a  two-phase  mixture  is  the  porosity,  a,  i.e.,  the  ratio  of  volume  occupied  by  the 
gas  phase  to  the  total  volume.  Ishii  (Ref.  21)  introduces  several  averages  which 
are  required  in  the  present  analysis.  Gough  (Ref.  22)  introduces  a  general 
weighting  function  g(y-x,i-t)  which  reflects  the  influence  of  remote  points  (y , t ) 
on  the  average  value  at  (x,t).  By  definition,  the  Gough  average  gives 


i 

Oil  V,t 


g  (  x ,t)  dx  dt  =  | 


(l) 


The  porosity  is  defined  by 


a(x,t)=  J  g(y-x,  r-t)  dy  dr 
V 

gas 

(2) 

The  weighting  function,  g,  plays  a  role  similar  to  the  state  density  functions 
(M^  M2,  Mg)  introduced  by  Ishii  (Ref.  21,  p.  65).  The  basic  time  average 
introduced  by  Ishii  (Ref.  21,  p.  68)  is  denoted  by  a  single  overbar  (ip),  and 
this  is  equivalent  to  Gough's  (Ref.  2)  unnormalized  average.  The  phase  average 
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denoted  by  a  double  overbar  (ij;)  is  related  to  ip  by 


=  xlf  I  r 

vi/  =  —  =  —  /  g  (  y- x  ,  t  -  t )  y/(y,  t)  dy  dT 

a  a  V 

vnas 

(3) 


Eq.  (3)  defines  the  average  of  a  gas  property,  ,  since  the  integral  is  carried 


out  over  the  region  occupied  by  the  gas  phase,  V 


gas 


In  Ishii’s  approach  the 


equivalent  average  is  obtained  by  integrating  only  over  the  time  interval  for 

which  the  gas  phase  is  present  at  the  space  point  x.  Finally,  the  mass  weighted 

th 

average  for  a  property  of  the  k  -phase  ^  is  defined  by 


tf  fvh  fiA 

^  '  pk  '  pk  (4) 


This  average  is  also  known  as  the  Favrd  average,  hence  the  superscript  F  is 

used.  This  is  a  very  convenient  average  to  use  in  turbulent  flow  since  density 

fluctuations  may  be  eliminated  formally.  It  should  be  noted  that  the  quantity 

_  _ 

is  the  partial  density  of  k  -phase  while  is  the  actual  density,  so  that 

the  mixture  density  is  given  by 


2 

-  I 

K=l 


'k^k 


(5) 


where  =  a,  and  =  1-a. 

In  the  following  equations,  the  Favrd  average  is  introduced  where  it  is 
appropriate,  and  phase  average  values  are  used  otherwise.  The  Favre-averaged 
velocity  vector  is  written  as 


=i  F 
U 


u 


(6) 
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and  on  all  other  variables  (e.g.,  e,  h,  etc.)  the  superscript  F  is  dropped  for 
convenience.  The  fluctuating  component  of  any  variable  is  denoted  with  a 
superscript  prime,  .  All  quantities  pertaining  to  the  solid  phase  are  denoted 

by  the  subscript  p.  The  derivation  of  the  equations  is  discussed  in  some  detail 
in  Ref.  (1)  and  will  not  be  repeated  here.  The  resulting  equations  are 
Gas  Phase  Continuity 


d(ap ) 
~d\ 


+  V  *  (  a/>  u )  <=  r(  +  M  j  g 


(7) 


Solid  Phase  Continuity 


d(l-a) 

~ at 


+  V  •  [  ( I  -  a  )  up] 


(8) 


where  the  mass  source,  T,  is  due  to  propellant  burning  and  M 
mass  addition  due  to  primer  discharge. 

Gas  Phase  Momentum 


ig 


is 


the  rate  of 


d  (ap  0) 

at 


+  V  •  (  a/>  UU  )  =  -  aVp  +  V-  ^q(tt  +  irT  )] 


-(l-a)^-  <F>  +  upr  +  M  jgUjg 


(9) 


Solid  Phase  Momentum 


4(|-a)PpgP] 

at 

+  V • [(I  -  a 


+  V  •  [  (  I  -a)/3pUpU 

n  Sp 

)R  +  (l-a)-rf- 

P 


p]  =  -  (  I  -  a)  Vp 

<F>  - upr 


(10) 
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In  the  above  equations,  tt  and  tt  are  the  average  stress  tensor  and  the 

turbulent  stress  tensor  in  the  gas  phase,  respectively,  1R  is  the  average 

T 

granular  stress  tensor,  and  tt  is  the  solid  phase  turbulent  stress  tensor. 

T  .  P 

For  the  present  time  tt^  will  be  neglected  because  there  is  insufficient 

information  available  to  construct  a  constitutive  relation  for  it.  The 
i 

quantity  <F>  is  the  interfacial  drag  force  and  is  the  primer  gas  velocity. 


Gas  Phase  Energy  Equation 

In  the  present  formulation  it  is  desirable  to  write  the  energy  equation 
in  terms  of  the  mass-averaged  static  enthalpy  h  because  of  simplifications 
in  the  turbulence  time-averaging. 


a(a  /?h) 


at 


+  V  •  (  a p U h )  =  -  V  •  [  a  (  q  +  q  T  )] 

( 


+  -  (ap)  +a<P+  ape  +  A  +  Mig  (  h,g  +  y  Uig-Uig 


(11) 


where  $  is  the  mean  flow  dissipation  defined  in  Eq .  (28),  e  is  the  turbulence 

kinetic  energy  dissipation  rate,  Eq.  (50),  tu  is  the  enthalpy  of  the  primer 
gas,  and  A  is  the  energy  transfer  term  between  the  solid  and  gas  phases, 

Eq.  (29). 

Turbulence  Kinetic  Energy  Equation 

In  the  present  work  a  turbulence  kinetic  energy  equation  and  a  specified 
length  scale  equation  has  been  utilized.  Following  the  derivations  of  Refs. 
(25-26)  and  Ref.  (1),  one  obtains 


+  a 


d(apk)  „ 

— t- +V-(ay5uk)  =  V  •  (  a  — —  V  k  ) 

at  crk 

^t[2d:d  -  ~(V-u)2]-  ~  pkV-U  -  pc  +  S 


(12) 
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where  is  an  interphase  turbulence  production  term,  Eq.  (68),  and  k  is  the 
mean  gas  phase  turbulence  kinetic  energy 


(13) 


Gas  Phase  Mixture  Molecular  Weight  and  Specific  Heat  Equations 

In  the  present  two-phase  flow  analysis  the  gas  phase  species  and  gasified 
propellant  species  mass  fractions  are  not  required.  Therefore,  in  order  to  limit 
computer  requirements  the  individual  species  mass  conservation  equations  are  not 
solved,  but  rather  only  total  gas  and  solid  phase  continuity  equations  are  solved 
Therefore,  it  is  necessary  to  consider  transport  equations  for  the  inverse 
mixture  molecular  weight  (Z)  and  the  specific  heat  at  constant  pressure  (c  )  : 


—  (a^Z)  +V-(ay5UZ)  =  V-[armVz]  +  ZpT 


(14) 


where  r  is  the  turbulent  exchange  coefficient  for  species  diffusion  which  is 
defined  from  a  knowledge  of  the  Schmidt  number  in  the  turbulent  flow  of  gas 
mixtures , 


Meff 

Sceff 


(15) 


and  is  generally  taken  as  a  constant,  Sc^^  =  0.9.  The  effective  viscosity 

is  defined  by  Eq .  (56).  Further,  is  the  inverse  molecular  weight  of 
the  propellant  and  T  is  the  mass  source  due  to  propellant  burning,  Eq .  (18). 

A  similar  transport  equation  may  be  derived  for  the  specific  heat  by 
assuming  that  the  species  specific  heats  are  independent  of  temperature: 


diap c  ) 

at 


+  V  •  (ap  Ucp) 


V-[“rmVcp] 


+  (cp)pr 


(16) 


where  (c  )  is 
P  P 


the  specific  heat  at  constant  pressure  of 


the  propellant. 
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Particle  Radius  Equation 

The  average  particle  radius,  r^ ,  is  required  as  a  function  of  spatial 


location  and  time  for  the  constitutive  relations  specified  below.  The  appro¬ 
priate  equation  including  turbulent  diffusion  is 


(17) 


where  the  relation  for  T,  Eq.  (18),  has  been  incorporated  in  order  to  cast  the 
equation  for  the  average  particle  radius  r^  into  the  above  form. 

Interfacial  Mass  Transfer 

Following  Ref.  (22),  the  mass  source  due  to  propellant  burning  may  be 
written 


P 


(18) 


where  S  is  the  average  particle  surface  area,  V  is  the  average  particle 
•  i  P 

volume,  and  <d>  is  the  average  regression  rate  of  the  solid  phase,  Eq .  (79). 
Stress  Tensors 

The  gas  phase  stress  tensor  assuming  a  Newtonian  fluid  is 

^  -  2/id  -  (  H-  -  O  v-  un 

V  3  B/  (19) 

where  is  the  bulk  viscosity  coefficient  and  I)  is  the  total  deformation  tensor 
(or  rate  of  strain  tensor)  given  by  (Ref.  21,  p.  164) 


D  =  D  .  +  D  - 
b  i 


(20) 
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where  ID,  is  the  bulk  deformation  tensor, 
b 

IDb  =-^[(Vu)  +  (VU)T]  (21) 

and  ID.  is  the  interfacial  deformation  tensor  which  is  difficult  to  model 
l 

except  for  a  dispersed  flow  (Ref.  21,  p.  165),  hence  it  must  be  neglected  at 
present.  The  turbulent  flow  stress  tensor  in  the  gas  phase  is  modeled 
using  an  isotropic  eddy  viscosity  formulation,  i.e.. 


ttt  =  -^u'u'  =  2fj.TD +pk)n  (22) 

where  k  is  the  turbulence  kinetic  energy,  Eq.  (13).  The  turbulent  viscosity 
must  be  determined  using  a  suitable  turbulence  model .  The  solid  phase 
granular  stress  tensor,  R,  is  modeled  by  assuming  an  isotropic  normal 
stress ,  i.e., 


R  5  RpI 

hence  in  the  solid  phase  momentum  Eq .  (10), 


(23) 


V*  [  ( I  -  a)  r]  =  V  [( I  -  a )  Rp  ]  (24) 

Heat  Flux  Vectors 

=  -*T 

The  mean  heat  flux  vector  q  and  the  turbulent  heat  flux  vector  q  in  a 

two-phase  flow  may  be  written  as  (Ref.  21,  p.  165) 


Q  = 


r  _  Va 
-*[VT-—  <T, 


(25) 


and 
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_ t  t  t  r  -  V  Q  __  _  *1 

qT  =  -  kt  I  V  T - —  (  T •  -  T  ) 

L  a  '  J  (26) 

=  T 

where  k  is  the  mean  thermal  conductivity,  Eq.  (48),  k  is  a  turbulent  conduc¬ 
tivity,  Eq.  (57),  and  T\  is  the  mean  temperature  at  the  interface  between  the 
phases.  For  the  present  time  T\  will  be  taken  as  the  average  between  the  gas 
temperature  and  the  particle  surface  temperature,  i.e. 


h  -T(T+TPS> 


(27) 


and  T  will  be  determined  from  the  solid  phase  heat  conduction  model, 
ps 

Mean  Flow  Dissipation 

The  mean  flow  dissipation  term  appearing  in  the  energy  equation,  Eq .  (11), 
is  defined  as 


$  =  2/i.  ED  :  D 


(-§-,* -Kb)(V.U>! 


(28) 


Interfacial  Energy  Transfer 


Following  Gough  (Ref.  22),  it  can  be  shown  that  the  interfacial  energy 
transfer  term  in  Eq.  (11)  is 


A  =  -  p(u-up)  •  Va 


+  (  I  -  a)  — 1 -  (  U-Up)  -  <F>‘  +  q-Va 
VP 


-  ( I  -  a  )  <  q  >  +  r  [ 

P 


I  _  _  _  _ 

h  .  +  —  (  u  -  u,J  •  ( u -  u_ 

comb  2  P  P 


>] 


(29) 


where  <q>1  is  the  interfacial  average  heat  transfer  between  the  gas  and  solid 
phases,  Eq.  (78),  and  h^^  is  the  energy  released  (per  unit  mass)  due  to 
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combustion  of  the  solid  propellant. 

Further  details  on  the  derivation  of  the  interphase  transfer  terms  in 
Eqs .  (7-17)  may  be  found  in  Refs.  (21-22),  and  a  summary  is  given  in  Ref.  1. 

In  the  present  analysis  Eqs.  (7-17)  are  solved  in  conjunction  with  Eqs.  (18-29) 

and  the  constitutive  relations  for  <d>^,  <F>*,  <q>^,  P,  K^,  k  and  which 

are  given  in  a  subsequent  section. 

Solid  Phase  Heat  Conduction  Equation 

Since  the  solid  particle  surface  temperature  is  desired  to  determine 
ignition,  the  propellant  burning  rate,  and  the  rate  of  heat  transfer  between  the 
gas  and  solid  phase,  a  transient  heat  conduction  equation  must  be  solved. 

Gough  (Ref.  2)  and  Kuo,  et  al.,  (Ref.  4)  assume  that  the  penetration  depth  of 
a  thermal  wave  into  the  propellant  grains  is  small  compared  to  the  grain  dimen¬ 
sions.  Then  it  is  permissible  to  use  a  one-dimensional  approximation  (planar 
for  cord  propellant  or  spherical  for  spherical  propellant  grains)  to  obtain  the 
particle  surface  temperature.  Following  the  motion  of  a  given  particle  (Kuo, 
et  al.,  Ref. 4),  the  heat  conduction  equation  for  a  spherical  particle  is 

/dfp\  _  Op_  ^(^Tp) 

\dt  )~  r2  a?2  oo) 


where  =  Tp  (r;  x,t)  is  the  phase-averaged  temperature  within  a  representative 

particle,  r  is  radial  coordinate  within  the  particle,  ot  is  the  thermal  diffus- 

ivity  of  the  solid  particles  [a  -  k  /P  (c  )  ],  and  (d/dt)~  denotes  the 
J  P  P  P  P  P  r 

Lagrangian  time  derivative  at  constant  r  within  the  particle.  Since  the  surface 
of  a  representative  burning  particle  is  receding  in  time  it  is  desirable  to 
employ  the  following  time-dependent  transformation  for  the  particle  radial 
coordinate  r: 


r 


i  0  <  £  <  l 


(31) 


Then  Eq.  (30)  becomes, 


dTp\  /  C  drp)  aTp 

dl  /£  \  rp  dl  j  d C 


(32) 
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where  the  quantity 


dr, 


Rs  "  dt 


=  <d> 


(33) 


may  be  identified  as  the  average  surface  regression  rate  for  the  particle, 
The  initial  condition  for  Eq.  (32)  is 


R  >  0. 
s  “ 


VC. i-o)  -  t 


po 


(34) 


The  boundary  conditions  are 


(C=  o,n  -  o 


at  C,=  0 


(35) 


h.  fie. 


(;  =  l,t)  -  hc(t)[f-fps]+qrad  +  k  i(Rs,p)  at 


(36) 


where  is  the  net  incident  radiation  heat  flux  normal  to  the  particle 

RAD 

surface,  k  is  the  thermal  conductivity  of  the  solid  particles,  and  ^(R^p) 
P  s 

is  the  heat  feedback  from  the  flame  identified  by  Gough  (Ref.  2,  p.  57). 

Assuming  that  the  gas  is  nearly  in  radiative  equilibrium  so  that  the  gas 

emissivity  is  unity,  and  that  radiation  emitted  by  other  particles  does  not 

influence  the  particle  in  question,  we  obtain 


’rad 


/  =A 

cr  (  T 


TPS> 


(37) 


where  e  is  the  particle  emissivity.  Other  authors  (e.g..  Refs.  2-4)  have 
P 

cast  Eq.  (37)  into  a  heat  transfer  coefficient  form,  so  Eq .  (36)  becomes 
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(38) 


if-jf-  (  C  =  *  hlO[T-TpJ  +  kp^(Rs.P) 

where  the  total  heat  transfer  coefficient  is 


h.  ■  hc  +Str(T  +Tps)(=f2  +  Tps> 


(39) 


The  convective  heat  transfer  coefficient,  h^  ,  is  specified  via  a  constitutive 
relation,  Eq.  (77).  An  expression  for  (j>  (Rg » p)  has  been  presented  by  Gough 
(Ref.  2)  for  a  planar  geometry  under  the  assumption  that  the  flame  zone 
surrounding  the  burning  particle  remains  quasi-s teady ,  and  that  the  convection 
and  radiation  heat  transfer  terms  in  Eq.  (36)  are  zero.  It  then  follows  that 


4> 


R, 


(Tps-fpo> 


(40) 


where  a  is  the  thermal  diffusivity  of  the  particles  and  T  is  the  undis- 
P  F  po 

turbed  temperature  far  from  the  particle  surface.  In  the  context  of  spherical 
particles,  T  would  be  taken  as  the  temperature  at  the  center  of  the  particle. 
This  procedure  should  be  sufficiently  accurate  in  view  of  the  other  assumptions 
made  in  obtaining  Eq .  (40). 

In  the  ALPHA2  code,  Eqs .  (34-37,  40)  have  been  utilized.  Solution  of 
this  solid  phase  heat  conduction  model  requires  special  consideration  since  it 
is  in  Lagrangian  form,  whereas  all  other  differential  conservation  equations 
are  in  Eulerian  form.  The  method  to  be  employed  in  the  present  analysis  will 
be  described  in  the  section  on  Solution  Procedure. 
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Constitutive  Relations 


The  necessary  constitutive  relations  include  a  gas  phase  equation  of  state,  a 
caloric  equation  of  state,  a  turbulence  length  scale  distribution,  the  molecular 
viscosity  and  thermal  conductivity,  the  so-called  form  functions  for  the  surface 
area  and  volume  of  the  solid  particles,  an  intergranular  stress  relation,  inter¬ 
phase  drag  and  heat  transfer  relations,  and  a  burning  rate  correlation  for  the  solid 
phase  combustion.  In  the  following,  the  double  overbar  (  )  is  dropped  for  simplicity. 

Gas  Equation  of  State 

The  Nobel-Abel  equation  of  state  will  be  used  for  the  gas, 

pRT 

pCI -pr,)  .  — -  5  pZl 
YYm 

where  is  the  universal  gas  constant,  is  the  gas  molecular  weight,  and  q  is  the 
covolume  factor,  which  is  composition  dependent.  Following  Gough  (Ref.  2)  an  arith¬ 
metic  average  will  be  used  for  q  based  upon  the  propellant  properties. 

The  caloric  equation  of  state  is  taken  as 


e 


cvT 


(42) 


where  is  dependent  on  the  gas  composition  but  not  temperature.  The  static  en¬ 
thalpy  is  then 


,  P  ZT 

h  =  e  +  —  =  c  T  + 


IP 


The  specific  heat  at  constant  pressure  is 


ah 

cp  =  aTlp  cv  +  z 


so  that  Eq .  (43)  may  be  written  as 


h  =  CdT  +  'PP 


(43) 


(44) 


(43) 
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Molecular  Viscosity,  Bulk  Viscosity,  and  Thermal  Conductivity  of  Gas 

The  molecular  viscosity  for  the  gas  is  determined  from  Sutherland’s  law, 


JL.  .  /  J  \3/2  ys, 

H- 0  h0  /  T  +  S, 


(46) 


where  =  110°K  for  air. 

The  bulk  viscosity  for  the  gas  will  be  assumed  to  be  zero  at  present, 


kb-o 


(47) 


The  thermal  conductivity  may  be  determined  from  a  relation  similar  to 
Sutherland’s  law,  e.g., 


K 


K 


0 


(  T  f/2  _V_S* 

ItJ  T  +  s2 


.  (48) 


where  S^  =  194°K  for  air. 


Turbulence  Model  Relations 

The  turbulent  viscosity  introduced  in  Eq .  (22)  is  obtained  from  the  Prandtl- 
Kolmogorov  relation,  viz. 


/*T 


H- 


(49) 


and  the  dissipation  rate  is  given  by 


€ 


_  3/2 
3 14  k 


(50) 


where  the  turbulence  length  scale,  £,  must  be  specified  consistent  with  the  expected 
turbulence  structure  in  the  two-phase  flow.  Following  Ref.  27  the  constant  will 
be  taken  as 


CT  =1.0 

k 


(51) 


The  very  large  fluid  accelerations  experienced  in  the  interior  ballistics 
problem  require  the  consideration  of  both  forward  and  reverse  transition  of 
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the  turbulent  flow  near  solid  surfaces.  There  are  two  options  available  for  model¬ 
ing  the  turbulence  near  a  wall.  In  the  first,  grid  point  resolution  normal  to  the 
wall  must  be  sufficient  to  define  the  viscous  sublayer,  in  which  case  the  boundary 
conditions  are  relatively  straightforward.  However,  the  difficulty  with  this  ap¬ 
proach  is  that  the  physics  of  low  Reynolds  number  (transitional)  turbulence  must  be 
modeled  in  a  reasonable  manner  by  the  governing  turbulence  equations  (e.g.,  Jones  and 
Launder,  Ref.  27).  An  alternative  approach  is  to  employ  a  less  refined  mesh  and 
force  the  turbulence  variables  to  yield  values  consistent  with  a  boundary  layer  wall 
function  formulation  at  the  first  grid  point  away  from  the  wall.  The  difficulty 
with  this  approach  is  that  the  validity  of  the  wall  function  formulation  is  question¬ 
able  under  the  rapidly  accelerating  transient  flow  situation  present  in  the  interior 
ballistics  problem.  Furthermore,  recent  experience  at  SRA  indicates  that  the  wall 
function  approach  may  be  inadequate  for  a  reacting  unsteady  flow  with  moving  coordi¬ 
nates.  Therefore,  a  transition  model,  which  was  sucessfully  used  by  Shamroth  and 
Gibeling  (Ref,  28)  in  a  time  dependent  airfoil  flow  field  analysis,  has  been  imple¬ 
mented  in  the  computer  code.  The  analysis  of  Ref.  28  follows  the  integral  turbulence 
energy  procedure  of  Refs.  29-31,  by  utilizing  a  turbulence  function  a^  ,  where 

0(=C//'2/2  (52) 

and  is  taken  as  a  function  of  a  turbulence  Reynolds  number  of  the  form 


— 1 

X) 

_ j 

/ 

f  (  R  T  )  1 

°,  ■  °0 

100 

/ 

1 . 0  +  6.66  o 

0 

- —  -  1 

100 

(53) 


where 


a0  =  .0115 

f(  Rr)  -  100.  Rr°  22  Rt  <  | 
f(RT)  =  68.1  Rt  +  614.3  Rt>40 


(54) 


and  a  cubic  curve  was  fit  for  values  of  R  between  1  and  40.  As  previously  discus¬ 
sed,  Ref.  29-31  utilized  an  integral  form  of  the  turbulence  kinetic  energy  and 
therefore  R^  was  defined  as  an  average  value. 
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(55) 


x/dy# 


In  the  present  effort  was  defined  as  the  local  ratio  of  turbulent  to  laminar  vis¬ 

cosity,  a  was  evaluated  via  Eq.  53  and  C  related  to  a,  via  Eq .  52. 
l  y  1 

The  effective  viscosity  is  defined  as  the  sum  of  the  laminar  and  turbulent 
viscosities 

Meff  =  M  +  (56) 


The  effective  conductivity  will  be  modeled  using  an  effective  Prandtl  number  obtained 
from  knowledge  of  turbulent  flows  of  gases  and  gas  mixtures,  i.e., 


eff 


=  K  +  K 


p^eff 


Pr 


eff 


and  Pre^f  =0.9  for  air. 


(57) 


Turbulence  Length  Scale 

For  the  evaluation  phase  of  the  present  effort,  the  turbulence  length  scale 
would  be  chosen  based  upon  known  steady  state  relations.  In  particular,  the  length 
scale  would  be  taken  as  the  minimum  of  the  length  scales  based  upon  the  local  average 
distance  between  solid  particles,  the  local  value  computed  from  turbulent  pipe  flow 
correlations,  and  that  from  turbulent  boundary  layer  length  scale  distributions  when 
close  to  the  wall. 

The  turbulent  pipe  flow  mixing  length  model  is  based  on  the  correlation  of 
Williamson  (Ref.  32), 


F( 


(58) 


where  r^  is  the  local  pipe  radius,  y  is  the  distance  from  the  point  in  question  to 

the  nearest  wall,  and  the  constant  C  is  taken  as 

w 

C  =  0.14  (59) 

w 


23 


The  near  wall  region  mixing  length  is  obtained  from  the  Van  Driest  model  (Ref.  33), 


£w  =  *y  ['  "  exp  (-y+/A+)] 


(60) 


where  k  is  the  von  Karman  constant  and  A  is  the  van  Driest  damping  coefficient, 


K  -  0.4 
A+  =  26.0 


(61) 


+ 


The  nondimensional  distance  y  is  defined  as 


>•  •  > « 


(62) 


and  the  friction  velocity  in  the  present  analysis  is  taken  as 


■  ft) 


1/2 


(63) 


where  the  local  shear  stress  is  obtained  from 


=  (  2D’ ID) 


1/2 


(64) 


The  average  center  to  center  distance  between  solid  spherical  particles 
assuming  dense  packing  is  approximately 


lr 


VP 

I  -a 


1/3 


(65) 


Therefore,  the  minimum  distance  between  particle  surfaces  is 


J?P  =  £>  Q  2fp 


(66) 
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and  the  turbulence  length  scale  will  be  assumed  proportional  to  the  distance  £ 


Ip  =  (3p£.p 


(67) 


The  constant  should  be  selected  based  upon  empirical  observation  of  the  turbulence 
structure  in  two-phase  flows  with  solid  particles. 

Interphase  Turbulence  Production 

The  interphase  turbulence  production  term  S  in  Eq .  (12)  may  be  inter- 

preted  as  the  production  of  turbulence  kinetic  energy  in  the  gas  phase  due  to 
gasification  of  solid  particles  (Ref.  1), 

Sk  ~  *psr  (68) 

where  kpg  represents  an  average  value  for  1/2  (u ' •  u')  at  the  gas-solid  inter¬ 
face.  It  is  not  known  how  to  specify  k^^  3t  the  present  time,  and  consequently 

until  further  information  becomes  available  k  =0 

ps 

Form  Functions 

The  surface  area  and  volume  of  particles  have  been  presented  by  Gough 
(Ref.  2)  for  a  variety  of  propellant  types.  In  the  present  work,  spherical 
propellant  grains  will  be  considered,  hence 


(69) 


where  r^  is  the  mean  particle  radius  at  a  given  point  in  space  and  time. 

Other  propellant  types  could  easily  be  considered  within  the  present  framework. 
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Intergranular  Stress  Relation 

A  stress  relation  for  granular  propellant  has  been  given  by  Gough  (Ref.  2), 
Koo  and  Kuo  (Ref.  3),  and  Kuo,  et  al. ,  (Ref.  4)  for  the  case  when  the  average 
stress  Rp  is  independent  of  the  loading  history: 


ac-a 


P?QP  (I -a)  a 


if 


a  <  ac 


R. 


0 


if 


a  >  a 


(70) 


where  is  a  critical  porosity  above  which  there  is  no  direct  contact  between 
particles,  and  a^  is  the  speed  of  sound  in  the  solid  phase  specified  on  input. 


Interphase  Drag  Relation 

-y  1 

The  average  steady  state  interphase  drag  <F>  appearing  in  the  momentum 
equations,  Eq .  (9-10),  was  obtained  from  correlations  for  nonfluidized  (packed) 
regions  and  fluidized  (dispersed)  regions.  Gough  (Ref.  2,  p.  48)  recommended 
the  following  drag  correlation,  valid  for  high  particle  Reynold's  numbers, 

Re  >>1,  which  is  based  on  Ergun!s  (Ref.  34)  correlation  for  a  packed  bed  and 
modified  by  Anderssen's  correlation  (Ref.  35)  as  the  bed  becomes  dispersed. 
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where  is  the  settling  porosity  and  is  given  by 


a.  = 


I  +  0.01986 


I  -  a. 


a. 


-I 


(73) 
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Also,  is  the  relative  velocity  between  the  gas  and  solid  phases,  i.e., 


V  U‘UP 


and  the  particle  Reynold’s  number  is 


(74) 


Re 


P 


(75) 


Gough  (Ref.  8)  has  noted  the  potential  source  of  error  due  to  the  presence  of 
relatively  few  propellant  grains  in  the  radial  direction  of  a  typical  medium 
caliber  weapon  (15-20  cm.  diameter) ,  since  the  available  correlations  are  based 
on  flows  through  long  columns  where  end  effects  are  negligible. 


Interphase  Heat  Transfer  Relation 

For  convective  heat  transfer  between  the  gas  and  solid  particles  in  interior 
ballistics  calculations,  numerous  correlations  have  been  recommended  (e.g., 

Refs.  2-4,  7).  Gough  (Ref.  2)  advocates  the  Gelperin-Einstein  correlation  (Ref.  36) 
for  the  interphase  heat  transfer  with  granular  propellant  in  both  fluidized  and 
nonfluidized  regions.  The  Nusselt  number  for  this  correlation  is  (Ref.  2) 


Nup  =  2.0  +  0.4  Rep/3  Prl/3 


(76) 


where  Pr  =  \ic^/k  and  k  is  the  gas  phase  thermal  conductivity,  Eq .  (48).  The 
heat  transfer  coefficient  in  Eq .  (36)  is  then 


(77) 


This  relation  is  considerably  simpler  than  the  Denton  and  Rowe-Claxton  correlations 
utilized  by  Kuo,  et  al.  (e.g.,  Refs.  3,  4),  and  may  yield  equally  reliable  predic¬ 
tions  in  view  of  the  large  variations  between  experimental  data  and  the  existing 
correlations  (Ref.  2). 

Finally,  the  interphase  heat  transfer  relation  required  in  the  energy  equa¬ 
tion  source  term,  Eq.  (29),  is 

<  q  >  h t  (  T  -  T ps)  (78) 

where  h  is  given  by  Eq.  (39). 
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Burning  Rate  Correlation 

The  steady  state  surface  regression  rate  (d  >  0)  is  given  by  (e.g.,  Ref.  2) 


d 


s 


(79) 


where  B^,  B2  and  n  have  known  constant  values.  The  phenomenon  of  erosive  burning 
is  assumed  to  be  an  acceleration  of  the  burning  rate  due  to  the  influence  of  con- 
vective  heat  transfer  on  the  heat  transfer  in  the  flame  zone.  The  Lenoir-Robillard 
(Ref.  37)  regression  rate  expression  is  utilized  for  this  effect, 


=  +  KEhcexp 


J 


(80) 


where  K^,  and  are  erosive  burning  constants,  determined  experimentally.  The  con¬ 
vective  heat  transfer  coefficient  is  then  obtained  from  Eqs.  (77,  78).  At  present 
only  the  steady  state  burning  relation,  Eq.  (79),  has  been  incorporated  into  the 
computer  code. 


Filler  Element  and  Projectile  Motion 

In  the  present  analysis  filler  elements  and  the  projectile  are  treated  dis¬ 
tinctly.  No  transverse  deformation  of  the  filler  elements  is  permitted  and  elements 
are  assumed  to  remain  planar;  therefore,  a  quasi-one-dimensional  lumped  parameter 
formulation  (e.g.,  Ref.  2)  may  be  employed  for  the  filler  elements.  The  appropriate 
equations,  which  have  been  stated  by  Gough  (Ref.  2)  are  repeated  here  for  complete¬ 
ness.  It  is  assumed  that  there  are  N  filler  elements  between  propellant  bed  and 
the  base  of  the  projectile,  with  the  projectile  denoted  as  element  (N+l) .  The 
required  properties  for  each  element  are  the  mass  (M^) ,  the  resistance  force  op¬ 
posing  motion  (F_^)  ,  an  internal  stress  (a_^)  ,  and  a  normal  wall  reaction  force  (Fw^) 
for  incompressible  elements  in  a  variable  area  tube.  The  cross-sectional  area  of 
each  element  is  assumed  to  be  equal  to  the  local  tube  area,  and  the  stress  in  an  in¬ 
compressible  element  is  assumed  to  be  isotropic. 

A  momentum  equation  is  then  written  for  one-half  of  element  i  together  with  one- 
half  of  element  (i-1)  in  order  to  describe  the  motion  of  the  interface  location,  z_^. 
There  results 


w 


—  M,  z(  =  A  j  o',  -  A0cr0- 


(81) 
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-|-,M|_I  +  Mi)z|  =  Aj  O',  -  A,  _ ,  cr  j_( 


T(Fi  +  F'-i> 


2  ‘  S  +  Fw,_,  ’ 


for  2  <  i  <  N 


(82) 


V  +  T-)  V,  '  '  '  (  F„»,  +  T-)  -  IT 


(83) 


In  this  section  the  stress  o_^  is  taken  as  positive  in  tension  following  Gough  (Ref. 

2),  and  the  tern  (-A  a  )  in  Eq.  (81)  is  the  force  exerted  on  the  first  filler  ele- 

o  o 

raent  by  the  gas  and  propellant  particles.  The  mass  of  the  projectile  is 

assumed  to  be  corrected  for  rotational  inertia;  if  I,  D  and  0  are  the  polar  moment 

D 

of  inertia,  the  tube  diameter  and  the  angle  of  rifling,  respectively,  it  follows 
that 

41  tan2# 

(84) 
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The  normal  wall  reaction  is  given  by 


0 


if  element  i  is  not  incompressible 


/  dA  \ 

(z,  -  Z  )  cr .  - )  if  element  i  is  incompressible 

1  +  1  1  1  \  d  Z  / . 


(85) 


Constitutive  data  must  be  provided  for  the  stress  a for  elastic  elements  or  for 

plastic  elements  in  a  state  of  loading  (i.e.,  ;  however,  for  rigid  elements 

or  plastic  elements  in  a  state  of  unloading  (z^z  ) ,  one  has 

i  1+ 1 


CTi  *°  ■  V 


(86) 


Finally,  for  an  incompressible  element,  i,  one  has  the  continuity  relation, 


V-,  =  ViVi 


(87) 
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Primer  Discharge 

The  primer  discharge  model  is  implemented  as  an  external  source  of  mass, 
momentum  and  energy  into  the  gas  phase  in  order  to  ignite  the  propellant  grains. 
The  general  functional  form  of  these  source  terms  is 

M  .  =  M  (  x .  t ) 

ig  igl*i 

%  =VV>  (88 

hig  =  hig<V<> 

In  practice  the  primer  source  is  added  only  in  a  small  region  near  either  the  tube 
centerline  for  a  center  core  igniter  or  the  breech  end  for  a  base  igniter. 
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COORDINATE  TRANSFORMATION 


The  set  of  governing  partial  differential  equations  which  model  the  physical 
processes  occurring  in  interior  ballistics  problems  was  presented  in  the  previous 
section.  For  generality  these  equations  were  written  in  vector  notation;  however, 
before  these  equations  can  be  incorporated  into  a  computer  code,  a  coordinate 
system  must  be  chosen.  The  governing  equations  can  then  be  cast  ip  a  form  reflect¬ 
ing  the  choice  of  the  coordinate  system.  The  coordinate  system  for  interior  bal¬ 
listics  calculations  must  have  the  ability  to  enlarge  the  physical  extent  of  the 
computational  domain  as  the  projectile  moves  through  the  gun  barrel,  and  it  must  be 
capable  of  treating  breech  geometries  in  which  the  tube  radius  is  a  function  of 
axial  distance.  Also,  a  time-dependent  transformation  is  required  to  adapt  the 
computational  mesh  to  regions  of  large  radial  gradients  appearing  during  the  igni¬ 
tion  phase  of  the  firing  cycle.  The  form  of  the  governing  equations  appearing  in 
Ref.  1  produces  mass  sources  due  to  time  truncation  error  when  radial  mesh  motion 
is  included  in  the  calculation.  This  observation  is  consistent  with  the  findings 
of  Thomas  and  Lombard  (Ref.  15),  and  appears  to  be  due  to  an  inconsistent  calcula¬ 
tion  of  the  overall  transformation  Jacobian  in  the  nonconservation  form  of  equations 
(Ref.  1). 

In  order  to  permit  consistent  calculation  of  the  Jacobian  in  a  moving  coordi¬ 
nate  system,  the  governing  equations  should  be  transformed  with  a  Jacobian  trans¬ 
formation  of  the  form 

yl  =  y'(x„x2,xs,t)  (89) 

T  =  t 

where  (x  ,  X3^  =  z ^  are  t^ie  original  cylindrical  polar  coordinates 

suitable  for  a  gun  tube.  The  velocity  components  remain  the  components ,  (bT  ?  U0 ,  U  ) 

in  the  (x^,  x^ ,  x^)  coordinate  directions  respectively.  The  new  independent  variables 
y~*  are  the  computational  coordinates  in  the  transformed  system.  The  coordinate 
system  requirements  for  interior  ballistics  applications  may  be  represented  by  a 
subset  of  the  general  transformation,  Eq .  (89), 


y'  =  y 1 ( x | ,  x3  5 1 ) 

(90a) 

y 2  =  y2  (x2  ) 

(90b) 

y3  =  y3(L,  *3, t) 

(90c) 
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which  is  a  general  axisymmetric  time-dependent  transformation.  For  interior  bal- 

2 

listics  flows  which  are  axisymmetric,  Eq.  (90b)  reduces  to  y  -  and  all  deriva- 

o 

tives  8/9y  are  assumed  to  be  zero.  The  transformation  (90)  with  the  axisymmetric 
flow  assumption  has  been  utilized  in  the  ALPHA2  computer  code. 

Application  of  the  Jacobian  transformation  requires  expansion  of  the  temporal 
and  spacial  derivatives  using  the  chain  rule,  i.e., 


=  J4.  +  y  vj  Hl  (91) 
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The  relations  Eq.  (91-93)  are  first  substituted  into  the  governing  equations  (7-17) 
written  in  cylindrical  polar  coordinates.  Then  the  resulting  equations  are  multi¬ 
plied  by  the  Jacobian  determinant  of  the  inverse  transformation. 
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and  the  equations  are  cast 
ing  relations. 


into  a  "semi-strong" 
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ay  * 
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conservation 


form  using 
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follow- 


(95) 
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and 


_aj_  +  £  ay],t 

dr  j=i  dyi 


=  0 


(96) 


where 


yj,i  -  JyV> 


a  i  __  .  -I 

y  ,t  =  Jy  ,t 


(97) 


The  semi-strong  conservation  form  implies  that  all  factors  involving  the  radial 
coordinate  r  =  x^  remain  as  they  were  before  the  Jacobian  transformation.  The 
resulting  equations  are  presented  in  Appendix  A. 

The  geometric  relations  Eq.  (95-96)  may  be  obtained  from  the  transformation 

/V  ♦  A  • 

i  i 

relations  for  y,^  and  y,  ^  in  terms  of  the  inverse  transformation  derivatives  (e.g., 
Ref.  15), 
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and 


*  j 

y ,» 


axk 


dr 


(99) 


The  projectile  motion  is  assumed  to  be  in  the  x^-direction ,  hence  the  follow- 

3  3  J 

ing  normalized  nonuniform  coordinate,  ri  (y  )  ,  is  introduced: 


y(y5)  = 


 1  ~~  z0 


Z I  Z0 


(100) 


where  z  is  the  Cartesian  location  of  the  breech  end  of  the  gun  barrel  and  z^  is  the 
o  3  i 

Cartesian  location  of  the  first  filler  element  and  y  is  an  equally  spaced  com- 

3 

putational  coordinate  having  a  value  from  1.0  to  y  .A  similar  transformation 

_  max 

can  be  implemented  for  the  radial  coordinate,  x^,  to  permit  the  tube  radius  to 
vary  with  axial  distance,  and  to  allow  a  concentration  of  grid  points  in  the  x^- 
direction  as  a  function  of  x^  and  time  t.  The  latter  feature  would  permit  the 
concentration  of  x^-direction  grid  points  in  the  manner  shown  in  Fig.  2,  to  account 
for  a  variation  of  the  boundary  layer  thickness  as  a  function  of  x^  and  t  as  the 
projectile  moves  through  the  barrel.  As  can  be  seen  in  Fig.  2  the  resulting 
coordinate  system  is  nonorthogonal ,  however,  such  a  system  is  already  encompassed 
by  the  general  transformation,  Eq .  (92).  Introducing  a  normalized  nonuniform 

coordinate , 


V  (y1) 


r  -  r 

ri~r 


o. 

0 


(101) 


where  r,  =  r.,  (xn)  is  the  radial  location  of  the  tube  wall,  r  =  r  (x0)  is  the 
113  i  o  o  3 

radial  location  of  a  centerbody  if  one  exists,  and  y  is  an  equally  spaced  compu¬ 
tational  coordinate  having  a  value  from  1.0  to  .  Of  course,  at  the 

^  max  i  i  33 

centerline  of  a  tube  r  =0.  The  functional  forms  of  n  (y  )  and  of  q  (y  )  are 

o 

arbitrary  and  can  be  chosen  such  that  the  packing  of  grid  points  in  the  or  X3- 
direction  is  achieved  in  the  regions  where  the  largest  gradients  are  expected. 
Presently  the  computer  code  allows  for  the  concentration  of  grid  points  to  occur 
by  means  of  Levy*s  (Ref.  38)  generalization  of  the  Roberts1  transformation  (Ref.  39). 
The  grid  points  can  be  concentrated  at  the  boundaries  of  the  computational  domain  or 
the  grid  points  can  be  concentrated  around  some  interior  location.  The  transforma¬ 
tion  equation  used  for  this  purpose  is 


34 


(102) 


V '  =  v'o*  |sinh{c[-lQnh(E^F)-H  ]  +  p} 


where  r|d  is  the  value  of  about  which  the  concentration  of  grid  points  is  centered 
and  the  values  of  A,  C,  D,  E,  F,  G  and  H  are  controlled  by  the  input  parameters, 

no’  C2’  T1  and  T2*  The  derivat:i-on  of  the  relationships  between  A,  C,  D,  E,  F,  G 
and  H  and  the  input  parameters  is  quite  lengthy  and  hence  only  the  results  are 
presented  here,  viz., 
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The  above  is  presented  only  for  completeness;  the  important  thing  to  note  is  the 

effect  that  nJ ,  t„,  T.  and  x„  have  on  the  physical  grid  spacing.  The  effect  of  t 
o  z  1  z 

is  to  regulate  the  sinh  portion  of  the  transformation,  while  x  and  regulate 

the  tanh  portion  of  the  transformation.  Note  that  different  parameter  values  may 

be  specified  for  each  of  the  coordinate  directions  x  in  order  to  control  the  grid 

spacing  in  each  of  those  directions;  x^  controls  the  physical  grid  spacing  at  the 

lower  end  (rp  .  )  of  the  computational  domain  while  x9  controls  the  spacing  at  the 
.  nun  z 


upper  end  (r^  ) 

max 


The  values  of  x^  and  are  subject  to  the  following  limitations 


-  I  <  TJ  <  0 
0  <  r2  <  I 
l2>0 


(113) 


In  order  to  see  how  the  input  parameters  effect  the  grid  spacing  it  is  instructive 
to  first  negate  the  effect  of  the  sinh  by  setting  *  0  and  to  investigate  the  ef¬ 
fect  that  x^  and  x^  have  on  the  transformation.  If  x^  =  0  and  T2  >  ^  packing 

will  occur  at  n3  (the  larger  the  value  of  x9  the  greater  the  packing)  while  if 
max  z. 

X  <  0  and  t2  =  0  packing  occurs  at  nJmin  (the  larger  the  value  of  |  t1  |  the  greater 
the  packing).  Zero  values  of  x  and  x  result  in  equal  grid  spacing  while  nonzero 
values  of  both  x^  and  x  result  in  packing  at  both  Trmin  and  n  T1=T2  and 

t 2~0  the  transformation  Eq.  (102)  is  equivalent  to  the  original  transformation  of 

Roberts  (Ref.  39),  and  the  parameter  X2  is  related  to  Roberts1  nonnalized  boundary 

1/2 

layer  thickness  parameter  a  by  x^  =  (1-a)  .  On  the  other  hand  if  the  effect  of 

the  tanh  is  negated  by  setting  both  x  and  x^  equal  to  zero,  the  effect  is  to  con- 

ct 

centrate  the  grid  points  about  n  only.  The  larger  the  value  of  t9  the  greater  the 

concentration.  Nonzero  values  of  t^,  and  x9  result  in  a  combination  of  the 

effects  of  the  sinh  and  the  tanh  transformations. 
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SOLUTION  PROCEDURE 


The  development  of  the  ALPHA2  computer  code  is  based  upon  an  axisymmetric 
version  of  the  highly  efficient  consistently  split,  linearized  block- implicit  solu¬ 
tion  procedure  (MINT)  for  the  compressible  Navier-Stokes  equations  developed  by 
Briley  and  McDonald  (Ref,  16-18),  and  subsequently  extended  to  multi-component, 
chemically  reacting,  turbulent  flows  by  Gibeling,  McDonald  and  Briley  (Ref.  19). 

This  procedure  solves  the  Navier-Stokes  equations  written  in  primitive  variables; 
in  the  MINT  procedure,  the  governing  equations  are  replaced  by  either  a  Crank- 
Nicholson  or  a  backward  time  difference  approximation.  Terms  involving  nonline¬ 
arities  at  the  implicit  time  level  are  linearized  by  Taylor  series  expansion  about 
the  known  time  level,  and  spatial  difference  approximations  are  introduced.  The 
result  is  a  system  of  two-dimensional  coupled  linear  difference  equations  for  the 
dependent  variables  at  the  unknown  or  implicit  time  level.  These  equations  are 
solved  by  the  Douglas-Gunn  (Ref.  40)  procedure  for  generating  ADI  schemes  as  per¬ 
turbations  to  fundamental  implicit  difference  schemes.  This  technique  leads  to 
systems  of  one-dimensional  coupled  linear  difference  equations  which  are  solved  by 
standard  block-elimination  methods,  with  no  iteration  required  to  compute  the  solu¬ 
tion  for  a  given  time  step.  An  artificial  dissipation  term  based  upon  either  a  cell 
Reynolds  number  criterion  or  the  rate  of  change  of  the  dependent  variable  may  be 
introduced  selectively  into  the  scheme  to  allow  calculations  to  be  performed  at 
high  local  values  of  the  cell  Reynolds  number. 

The  use  of  an  implicit  solution  procedure  requires  that  equation  coupling 
and  linearization  be  considered.  Both  of  these  questions  are  reviewed  in  detail 
by  McDonald  and  Briley  (Ref.  41)  and  Briley  and  McDonald  (Ref.  18).  These  authors 
have  argued  that  for  a  given  grid  the  errors  arising  from  time  linearization  of  the 
nonlinear  terms  at  the  unknown  time  level  should  be  no  greater  than  the  discretiza¬ 
tion  errors.  Also,  reduction  of  the  time  step  is  the  preferred  way  of  reducing  the 
linearization  error  since  transient  accuracy  is  thereby  improved.  Linearization  by 
Taylor  series  expansion  in  time  about  the  known  time  level  introduces  errors  no 
greater  than  those  due  to  the  differencing  (Refs.  41  and  18),  and  this  approach  has 
been  employed  in  the  ALPHA2  code.  The  formal  linearization  process  results  in  a 
system  of  coupled  equations  in  order  to  retain  second-order  temporal  accuracy. 

The  system  of  coupled  equations  at  the  implicit  time  level  is  solved  efficiently 
using  a  standard  block  elimination  matrix  inversion  scheme.  In  the  present  problem, 
the  strong  coupling  effects  among  the  governing  equations  dictate  the  use  of  the 
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block  coupled  equation  approach.  However,  weakly  coupled  equations  would  probably 
be  solved  in  a  decoupled  manner  in  order  to  reduce  computer  time  and  storage 
requirements . 

The  principal  partial  differential  equations  which  will  be  solved  via  the 
MINT  technique  are:  gas  and  solid  phase  continuity,  gas  and  solid  phase  momenta, 
gas  phase  energy,  gas  phase  turbulence  kinetic  energy,  gas  phase  mixture  molecular 
weight  and  specific  heat  equations  and  the  particle  radius  equation.  The  constitu¬ 
tive  relations  required  to  close  the  above  system  of  equations  have  been  specified 
above.  The  solid  phase  heat  conduction  equation  is  the  only  differential  equation 
which  requires  special  treatment  because  it  is  a  Lagrangian  equation. 

The  scheme  devised  for  solution  of  the  solid  phase  heat  conduction  equation 
is  unique  since  it  does  not  involve  the  use  of  marker  particles  introduced  by  other 
authors  (e.g. ,  Ref.  2).  This  is  possible  because  the  equation  is  a  simple  heat 
conduction  equation  for  a  representative  solid  particle  moving  at  a  velocity 
which  is  known  at  the  completion  of  a  given  time  step.  The  necessary  boundary 
conditions,  Eqs.  (35-38),  provide  information  about  the  environment  through  which 
the  particle  is  moving  in  the  form  of  a  heat  transfer  coefficient,  Eq .  (39).  The 

procedure  to  be  used  assumes  that  at  time  tn~^  the  representative  particle  has  moved 

.  .  ,  .  -Hi+1  ,  n+1  n+1  n+1  .  .  n  ,  .  . 

to  the  grid  point  x.  .  ,  =  (xn  x0  from  a  location  at  time  t  which  is 

1  >  3  ,  R  1  •  ,  3 .  5  3.  ) 

1  1  k 

~>n+ 1  -hi 

determined  from  the  known  absolute  particle  velocities  v  and  v  ,  i.e.,  if  s  is 

p  p  p 

the  particle  position  vector  relative  to  an  inertial  reference  frome,  we  have 

ds_ 

v_  (114) 


dt 


and  application  of  the  variable  time  differencing  scheme  yields. 


n+l  _ 

5p  -sp 


+[/3vp"*'+(,-/3)vpn]At 


(115) 


where  6=1  for  backward  time  differencing  and  6=1/2  for  Crank-Nicholson  (centered) 
time  differencing.  In  the  present  scheme,  is  assumed  to  be  the  grid  point  lo¬ 
cation  ,  and  Eq .  (115)  is  then  solved  for  sU.  Because  the  grid  is  moving  it  is 

i,J>k  p 

necessary  to  interpolate  to  find  the  required  value  of  the  particle  velocity  at 

n  .  -Hi+1  -Hi  -hi /->n+l  n.  , .  .  „  /oc\ 

time  t  at  space  point  x.  .  .  ,  v  =  v  (x.  .  .  ,  t  ).  The  boundary  condition,  Eq.  (3b), 

P  P  i,3  A 

may  then  be  written  as 
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(116) 


k  p  Tp 


rp  ac  (^l,tn*  )  */3[h,(tn*,)(f"*'-  fp"s4l)+  k^lR^'.p"*')] 
+  (l-/3)[ht(ln)(T"-Tpn5  )  +  kp^(Bsn,p")] 


n 


The  desired  properties  h  (t  )  ,  T  ,  T  ,  and  $ (R  ,  p  )  are  understood  to  be 

-kT  Ps  s 

evaluated  at  the  point  s^,  and  these  will  be  evaluated  by  interpolation  utilizing 

values  at  time  t  at  the  four  grid  points  surrounding  point  sR. 

P 

Finally,  the  governing  equation  (32)  and  boundary  conditions,  Eq .  (35)  and 

(116) ,  may  be  written  in  finite  difference  form.  The  resulting  tridiagonal  matrix 
is  easily  inverted  using  Gaussian  elimination  to  yield  the  temperature  distribution 
within  the  particle.  Another  approximate  solution  technique  could  be  incorporated 
at  a  later  time  in  order  to  reduce  the  computer  requirements  for  the  particle  heat 
conduction  model. 


Initial  and  Boundary  Conditions 

The  initial  conditions  for  the  first  phase  of  two-dimensional  calculations 
will  consist  of  a  description  of  the  fluidized  state  of  the  flow  in  a  gun  barrel 
after  ignition  is  complete  and  the  projectile  motion  has  begun.  Typically,  the 
necessary  data  would  be  produced  from  an  existing  one-dimensional  interior  bal¬ 
listics  computer  code,  and  would  then  be  extended  over  the  two-dimensional  computa¬ 
tional  domain  by  applying  a  correction  for  the  wall  boundary  layers.  Provisionally, 
the  boundary  layer  integral  method  adopted  by  Gough  (Ref.  2)  would  be  utilized  to 
determine  the  boundary  layer  thickness  and  velocity  profile. 

The  boundary  conditions  to  be  applied  would  be  no-slip  gas  velocities  on 
solid  surfaces  and  conventional  symmetry  conditions  at  the  tube  centerline.  The 
breech  would  be  assumed  to  be  stationary,  and,  of  course,  the  projectile  and  filler 
elements  would  be  allowed  to  move.  The  wall  pressure  would  be  determined  by  employ¬ 
ing  the  normal  gas  momentum  equation  written  at  the  wall.  The  surface  temperature 
would  be  determined  by  incorporating  a  barrel  heat  conduction  model  coupled  to  the 
gas  heat  transfer  at  the  wall.  For  simplicity,  heat  conduction  in  the  barrel  would 
be  assumed  to  be  primarily  in  the  radial  direction.  The  porosity  at  a  wall  would  be 
determined  from  either  the  solid  phase  continuity  equation,  Eq .  (8),  or  a  zero  normal 
derivative  condition.  The  solid  phase  velocities  at  a  wall  would  be  determined 
from  either  the  no-slip  condition  or  the  solid  phase  momentum  equations  written 
at  the  wall. 
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The  appropriate  boundary  condition  for  the  inverse  mixture  molecular  weight 
(Z) ,  specific  heat  (c^) >  and  particle  radius  (r  )  a  solid  wall  is  zero  normal 
derivative,  i.e., 


(117) 


This  follows  from  the  definition  of  these  quantities  as  mass  weighted  averages,  and 


the  assumption  that  the  individual  species  diffusion  velocities  normal  to  the  wall 
as  determined  from  Fick's  law  must  be  zero;  that  is,  (3m. /3n)  =  0  where  m  is  a 
species  mass  fraction. 

The  solid  particles  which  reach  the  wall  will  be  assumed  to  be  in  equilibrium 
with  the  gas  phase,  thus 


(118) 
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Figure  1.  Scht’in  it  ic  illustration  of  pun  breech  region  prior  to  firing 
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APPENDIX 


The  governing  conservation  equations  in  cylindrical-polar  coordinates  are 
transformed  using  the  Jacobian  transformation, 


yj  =  yj(X|,  x2,  x3,t) 

T  =  t 


(A-l) 


where  (x^,  x^,  x^)  =  (r,  G,  z) .  The  resulting  equations  may  be  written  in  the 
following  compact  form:. 


d{  JW  ) 
dr 


dyJ 
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(A- 3) 


Further,  the  coefficients  B^,  y^,  are  given  by 


~  r  >  @z  ~  r  »  ^3 =  1 

Y\  "  1  >  Yz  ~~  ~T  ’  Y*  =  1 


(A-4) 


and  m  -  1  for  all  equations  except  the  x„-  direction  momentum  equations  (gas  and 
solid  phase),  for  which  m  =  2.  The  quantity  a  is 


a  =  a  for  gas  phase  equations 

a  =  I  -  a  f°r  solid  phase  equations 


(A- 5) 
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The  vector  variables  used  in  Eq.  (A-2)  are  defined  as 


W  = 


where  n 
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Note  that  the  velocity  components  (U^,  U^,  U^)  and  (U  ,,  U  ,  U  «)  are  the 
cylindrical-polar  velocity  components,  and  ..  and  .  are  the  gas-phase  stress 
tensor  and  the  solid  phase  intergranular  stress  tensoif,  respectively,  written  in 
cylindrical-polar  coordinates.  The  molecular  and  turbulent  stress  tensors, 

Eqs .  (19)  and  (22),  may  be  written  as 


Tij  =  2/xeff  Dij 


2  M,„<v-U)8i|  + 4(KbV.U  -p k)SM 


3  '  eff  IJ  3  *  "B'  r-'~  jj 

and  the  rate  of  strain  tensor  components  in  cylindrical-polar  coordinates  are 
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(A-12) 
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and 


(A-13) 


Also,  an  isotropic  intergranular  stress  tensor  is  assumed, 


(A-14) 


where  the  constitutive  relation  for  R  is  given  by  Eq.  (70).  Therefore,  the  inter¬ 
granular  stress  becomes  a  normal  strels  as  is  the  gas  pressure.  The  derivatives 
required  in  Eqs.  (A-12,  A-13)  must  be  expressed  in  terms  of  the  computational 
coordinates  y3  using  the  chain  rule,  Eq.  (92). 

-y 

Finally,  the  vector  S  contains  source  terms  and  certain  dif f erential^terms 
which  do  not  conform  to  the  basic  structure  of  Eq.  (A-2) ,  and  the  vector  C  contains 
the  additional  curvature  terms  due  to  the  cylindrical-polar  coordinate  system. 
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where 


I  ..2  I  — 

r  a/°U2  r  a  r22 

-  a/?  U|  U2 
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0 

0 

0 


C  = 


0 

0 


0 

-f<'-a>P pup22  -  T11"0'^ 
--rl'-alPpUp.Upj 
0 


0 


(A-16) 


M 


S  ' 

-  (I  -  a)  -rp  <F>  +  Up  r 


(A-17) 


The  quantities  r,  <?>  1,  <d>  1,  and  A  are  defined  in  Eqs.  (18),  (71),  (79)  and 

(29),  respectively. 
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